#include "ManagedLasoProvider.h"

/* Table of constant values */

static doublereal c_b2 = 10.;
static integer c__1 = 1;
static integer c__0 = 0;
static doublereal c_b72 = 0.;


/* ------------------------------------------------------------------ */

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::diwla_(S_fp
   op, S_fp iovect, integer *n, integer *nband, doublereal *xl, doublereal *xr,
   integer *nfig, integer *nperm, doublereal *val, integer *nmvec, integer
   *maxvec, doublereal *vec, integer *nblock, integer *maxop, integer *maxj,
   integer *nop, doublereal *p1, doublereal *p0, doublereal *res, doublereal
   *tau, doublereal *otau, doublereal *t, doublereal *alp, doublereal *bet,
   doublereal *s, doublereal *p2, doublereal *bound, doublereal *atemp,
   doublereal *vtemp, doublereal *d__, integer *ind, logical *raritz,
   doublereal *deltal, doublereal *deltar, doublereal *eps, integer *ierr)
{
  /* System generated locals */
  integer vec_dim1, vec_offset, p0_dim1, p0_offset, p1_dim1, p1_offset, p2_dim1,
    p2_offset, t_dim1, t_offset, alp_dim1, alp_offset, bet_dim1, bet_offset,
    s_dim1, s_offset, i__1, i__2, i__3, i__4;
  doublereal d__1, d__2, d__3, d__4;

  /* Builtin functions */
  //double   pow_dd(doublereal *, doublereal *);

  /* Local variables */
  static logical done;
  /*extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
  integer *);*/
  static doublereal tola, temp, tolg, tmin, tmax, tarr[1];
  static integer numl;
  static logical test;
  static integer numr;
  static doublereal utol;
  /*extern doublereal dnrm2_(integer *, doublereal *, integer *);*/
  static integer i__, j, k, l, m;
  //extern /* Subroutine */ int dscal_(integer *, doublereal *, doublereal *, 
  // integer *);
  static integer ngood;
  //extern /* Subroutine */ int dmvpc_(integer *, doublereal *, integer *, 
  // integer *, doublereal *, integer *, doublereal *, doublereal *, 
  // doublereal *);
  static doublereal anorm;
  //extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *);
  static doublereal dzero[1];
  //extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
  // integer *, doublereal *, integer *);
  static integer i1;
  static doublereal epsrt, rnorm;
  static integer ii, ng;
  //extern /* Subroutine */ int dlaeig_(integer *, integer *, integer *, 
  // integer *, doublereal *, doublereal *, integer *, doublereal *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, doublereal *);
  static integer nl;
  //extern /* Subroutine */ int dlager_(integer *, integer *, integer *, 
  // doublereal *, doublereal *, doublereal *), dlaran_(integer *, 
  // doublereal *);
  static doublereal betmin, alpmin, betmax, alpmax;
  static integer ntheta, number, nstart;
  //extern /* Subroutine */ int dortqr_(integer *, integer *, integer *, 
  // doublereal *, doublereal *), dvsort_(integer *, doublereal *, 
  // doublereal *, integer *, doublereal *, integer *, integer *, 
  // doublereal *);
  static doublereal axl, axr;



  /* DIWLA IMPLEMENTS THE LANCZOS ALGORITHM WITH SELECTIVE */
  /* ORTHOGONALIZATION. */

  /*   NPERM   THE NUMBER OF PERMANENT VECTORS (THOSE EIGENVECTORS */
  /*     INPUT BY THE USER OR THOSE EIGENVECTORS SAVED WHEN THE */
  /*     ALGORITHM IS ITERATED).  PERMANENT VECTORS ARE ORTHOGONAL */
  /*     TO THE CURRENT KRYLOV SUBSPACE. */

  /*   NOP   THE NUMBER OF CALLS TO OP. */

  /*   P0, P1, AND P2   THE CURRENT BLOCKS OF LANCZOS VECTORS. */

  /*   RES   THE (APPROXIMATE) RESIDUAL NORMS OF THE PERMANENT VECTORS. */

  /*   TAU AND OTAU   USED TO MONITOR THE NEED FOR ORTHOGONALIZATION. */

  /*   T   THE BAND MATRIX. */

  /*   ALP   THE CURRENT DIAGONAL BLOCK. */

  /*   BET   THE CURRENT OFF DIAGONAL BLOCK. */

  /*   BOUND, ATEMP, VTEMP, D  TEMPORARY STORAGE USED BY DLAEIG. */

  /*   S   EIGENVECTORS OF T. */

  /*   RARITZ   RETURNED AS .TRUE. IF A FINAL RAYLEIGH-RITZ PROCEDURE */
  /*     IS TO BE DONE. */

  /*   DELTAL   RETURNED AS THE EXCLUDED EIGENVALUE CLOSEST TO XR. */
  /*     USED TO ESTIMATE THE ACCURACY OF THE COMPUTED EIGENVALUES. */

  /*   DELTAR   RETURNED AS THE EXCLUDED EIGENVALUE CLOSEST TO XR. */
  /*     USED TO ESTIMATE THE ACCURACY OF THE COMPUTED EIGENVALUES. */


  /* INTERNAL VARIABLES USED */


  /*   J   THE CURRENT DIMENSION OF T.  (THE DIMENSION OF THE CURRENT */
  /*     KRYLOV SUBSPACE) */

  /*   NGOOD   THE NUMBER OF GOOD RITZ VECTORS (GOOD VECTORS */
  /*     LIE IN THE CURRENT KRYLOV SUBSPACE). */

  /*   NUMBER = NPERM + NGOOD */

  /*   ANORM   AN EXTIMATE OF THE NORM OF THE MATRIX. */

  /*   EPS   THE RELATIVE MACHINE PRECISION. */

  /*   UTOL   THE USER TOLERANCE. */

  /*   DZERO IS AN ARRAY OF SIZE 1 WITH VALUE ZERO.  USED FOR CONSISTENCY OF */
  /*         VARIABLE TYPING IN CALLS TO DCOPY. */

  /*   TARR  IS AN ARRAY OF LENGTH 1 USED FOR CONSISTENCE OF VARIABLE TYPING */
  /*         IN CALLS TO DLAEIG. */

  /*  INITIALIZE */

  /* Parameter adjustments */
  p2_dim1 =  *n;
  p2_offset = 1+p2_dim1 * 1;
  p2 -= p2_offset;
  p0_dim1 =  *n;
  p0_offset = 1+p0_dim1 * 1;
  p0 -= p0_offset;
  p1_dim1 =  *n;
  p1_offset = 1+p1_dim1 * 1;
  p1 -= p1_offset;
  t_dim1 =  *nband;
  t_offset = 1+t_dim1 * 1;
  t -= t_offset;
  --val;
  vec_dim1 =  *nmvec;
  vec_offset = 1+vec_dim1 * 1;
  vec -= vec_offset;
  bet_dim1 =  *nblock;
  bet_offset = 1+bet_dim1 * 1;
  bet -= bet_offset;
  alp_dim1 =  *nblock;
  alp_offset = 1+alp_dim1 * 1;
  alp -= alp_offset;
  s_dim1 =  *maxj;
  s_offset = 1+s_dim1 * 1;
  s -= s_offset;
  --res;
  --tau;
  --otau;
  --bound;
  --atemp;
  --vtemp;
  --d__;
  --ind;

  /* Function Body */
  dzero[0] = 0.;
  rnorm = 0.;
  if (*nperm != 0)
  {
    /* Computing MAX */
    d__1 =  - val[1], d__2 = val[ *nperm];
    rnorm = max(d__1, d__2);
  }
  epsrt = sqrt(*eps);
  *nop = 0;
  number =  *nperm;
  *raritz = FALSE_;
  /* Computing MAX */
  d__3 = (doublereal)( - ((real)(*nfig)));
  d__1 = (doublereal)((real)(*n)) **eps, d__2 = pow_dd(&c_b2, &d__3);
  utol = max(d__1, d__2);
  /* Computing MAX */
  d__1 =  *xr, d__2 =  - (*xl);
  temp = utol * max(d__1, d__2);
  axl =  *xl + temp;
  /* Computing MAX */
  d__1 = axl, d__2 =  *xr - temp;
  axr = max(d__1, d__2);
  j =  *maxj;

  /* ------------------------------------------------------------------ */

  /* ANY ITERATION OF THE ALGORITHM BEGINS HERE. */

  L30: i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    temp = dnrm2_(n, &p1[i__ *p1_dim1 + 1], &c__1);
    if (temp == 0.)
    {
      dlaran_(n, &p1[i__ *p1_dim1 + 1]);
    }
    /* L50: */
  }
  if (*nperm == 0)
  {
    goto L70;
  }
  i__1 =  *nperm;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    tau[i__] = 1.;
    otau[i__] = 0.;
    /* L60: */
  }
  L70: i__1 =  *n **nblock;
  dcopy_(&i__1, dzero, &c__0, &p0[p0_offset], &c__1);
  i__1 =  *nblock **nblock;
  dcopy_(&i__1, dzero, &c__0, &bet[bet_offset], &c__1);
  i__1 = j **nband;
  dcopy_(&i__1, dzero, &c__0, &t[t_offset], &c__1);
  m =  *maxvec + 2;
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    dcopy_(&j, dzero, &c__0, &s[i__ *s_dim1 + 1], &c__1);
    /* L75: */
  }
  test = TRUE_;
  ngood = 0;
  betmax = 0.;
  tmin = 1e30;
  tmax =  - 1e30;
  numl = 1;
  numr = 1;
  j = 0;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION TAKES A SINGLE BLOCK LANCZOS STEP. */

  L80: j +=  *nblock;

  /* THIS IS THE SELECTIVE ORTHOGONALIZATION. */

  if (number == 0)
  {
    goto L110;
  }
  i__1 = number;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (tau[i__] < epsrt)
    {
      goto L100;
    }
    test = TRUE_;
    tau[i__] = 0.;
    if (otau[i__] != 0.)
    {
      otau[i__] = 1.;
    }
    i__2 =  *nblock;
    for (k = 1; k <= i__2; ++k)
    {
      temp =  - ddot_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p1[k *p1_dim1 + 1],
        &c__1);
      daxpy_(n, &temp, &vec[i__ *vec_dim1 + 1], &c__1, &p1[k *p1_dim1 + 1],
        &c__1);

      /* THIS CHECKS FOR TOO GREAT A LOSS OF ORTHOGONALITY BETWEEN A */
      /* NEW LANCZOS VECTOR AND A GOOD RITZ VECTOR.  THE ALGORITHM IS */
      /* TERMINATED IF TOO MUCH ORTHOGONALITY IS LOST. */

      if ((d__1 = temp * bet[k + k * bet_dim1], abs(d__1)) > (doublereal)((real)
        (*n)) *epsrt * anorm && i__ >  *nperm)
      {
        goto L300;
      }
      /* L90: */
    }
    L100: ;
  }

  /* IF NECESSARY, THIS REORTHONORMALIZES P1 AND UPDATES BET. */

  L110: if (!test)
  {
    goto L160;
  }
  test = FALSE_;
  dortqr_(n, n, nblock, &p1[p1_offset], &alp[alp_offset]);
  if (j ==  *nblock)
  {
    goto L160;
  }
  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (alp[i__ + i__ * alp_dim1] > 0.)
    {
      goto L130;
    }
    m = j - (*nblock << 1) + i__;
    l =  *nblock + 1;
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      bet[i__ + k * bet_dim1] =  - bet[i__ + k * bet_dim1];
      t[l + m * t_dim1] =  - t[l + m * t_dim1];
      ++m;
      --l;
      /* L120: */
    }
    L130: ;
  }

  /* THIS IS THE LANCZOS STEP. */

  L160: (*op)(n, nblock, &p1[p1_offset], &p2[p2_offset]);
  ++(*nop);
  (*iovect)(n, nblock, &p1[p1_offset], &j, &c__0);

  /* THIS COMPUTES P2=P2-P0*BET(TRANSPOSE) */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      d__1 =  - bet[i__ + k * bet_dim1];
      daxpy_(n, &d__1, &p0[k *p0_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1], &c__1)
        ;
      /* L170: */
    }
    /* L180: */
  }

  /* THIS COMPUTES ALP AND P2=P2-P1*ALP. */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      i1 = i__ - k + 1;
      alp[i1 + k * alp_dim1] = ddot_(n, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k
        *p2_dim1 + 1], &c__1);
      d__1 =  - alp[i1 + k * alp_dim1];
      daxpy_(n, &d__1, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1)
        ;
      if (k != i__)
      {
        d__1 =  - alp[i1 + k * alp_dim1];
        daxpy_(n, &d__1, &p1[k *p1_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1],
          &c__1);
      }
      /* L190: */
    }
    /* L200: */
  }

  /*  REORTHOGONALIZATION OF THE SECOND BLOCK */

  if (j !=  *nblock)
  {
    goto L220;
  }
  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      temp = ddot_(n, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1);
      d__1 =  - temp;
      daxpy_(n, &d__1, &p1[i__ *p1_dim1 + 1], &c__1, &p2[k *p2_dim1 + 1], &c__1)
        ;
      if (k != i__)
      {
        d__1 =  - temp;
        daxpy_(n, &d__1, &p1[k *p1_dim1 + 1], &c__1, &p2[i__ *p2_dim1 + 1],
          &c__1);
      }
      ii = i__ - k + 1;
      alp[ii + k * alp_dim1] += temp;
      /* L210: */
    }
    /* L215: */
  }

  /* THIS ORTHONORMALIZES THE NEXT BLOCK. */

  L220: dortqr_(n, n, nblock, &p2[p2_offset], &bet[bet_offset]);

  /* THIS STORES ALP AND BET IN T. */

  i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m = j -  *nblock + i__;
    i__2 =  *nblock;
    for (k = i__; k <= i__2; ++k)
    {
      l = k - i__ + 1;
      t[l + m * t_dim1] = alp[l + i__ * alp_dim1];
      /* L230: */
    }
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      l =  *nblock - i__ + k + 1;
      t[l + m * t_dim1] = bet[k + i__ * bet_dim1];
      /* L240: */
    }
    /* L250: */
  }

  /* THIS SHIFTS THE LANCZOS VECTORS */

  i__1 =  *nblock **n;
  dcopy_(&i__1, &p1[p1_offset], &c__1, &p0[p0_offset], &c__1);
  i__1 =  *nblock **n;
  dcopy_(&i__1, &p2[p2_offset], &c__1, &p1[p1_offset], &c__1);
  i__1 = j -  *nblock + 1;
  dlager_(&j, nband, &i__1, &t[t_offset], &tmin, &tmax);
  /* Computing MAX */
  d__1 = max(rnorm, tmax), d__2 =  - tmin;
  anorm = max(d__1, d__2);
  if (number == 0)
  {
    goto L265;
  }

  /* THIS COMPUTE THE EXTREME EIGENVALUES OF ALP */

  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, 
    &tmax);
  alpmin = tarr[0];
  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  dlaeig_(nblock, nblock, nblock, nblock, tarr, tarr, nblock, &p2[p2_offset],
    &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, &tmax);
  alpmax = tarr[0];

  /*  THIS COMPUTES ALP = BET(TRANSPOSE)*BET */

  L265: i__1 =  *nblock;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 = i__;
    for (k = 1; k <= i__2; ++k)
    {
      l = i__ - k + 1;
      i__3 =  *nblock - i__ + 1;
      alp[l + k * alp_dim1] = ddot_(&i__3, &bet[i__ + i__ * bet_dim1], nblock,
        &bet[k + i__ * bet_dim1], nblock);
      /* L270: */
    }
    /* L275: */
  }
  if (number == 0)
  {
    goto L290;
  }

  /*  THIS COMPUTES THE SMALLEST SINGULAR VALUE OF BET */

  dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  d__1 = anorm * anorm;
  dlaeig_(nblock, nblock, &c__1, &c__1, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &c_b72, 
    &d__1);
  betmin = sqrt(tarr[0]);

  /* THIS UPDATES TAU AND OTAU. */

  i__1 = number;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* Computing MAX */
    d__1 = alpmax - val[i__], d__2 = val[i__] - alpmin;
    temp = (tau[i__] *max(d__1, d__2) + otau[i__] *betmax +  *eps * anorm) /
      betmin;
    if (i__ <=  *nperm)
    {
      temp += res[i__] / betmin;
    }
    otau[i__] = tau[i__];
    tau[i__] = temp;
    /* L280: */
  }

  /* THIS COMPUTES THE LARGEST SINGULAR VALUE OF BET */

  L290: dcopy_(nblock, dzero, &c__0, &p2[p2_offset], &c__1);
  d__1 = anorm * anorm;
  dlaeig_(nblock, nblock, nblock, nblock, &alp[alp_offset], tarr, nblock,
    &p2[p2_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &c_b72,
    &d__1);
  betmax = sqrt(tarr[0]);
  if (j <=  *nblock << 1)
  {
    goto L80;
  }

  done = TRUE_;
  goto L310;

  /* THIS INDICATES UNEXPECTED LOSS OF ORTHOGONALITY. */

  L300: j -=  *nblock;
  *ierr =  - 8;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES SOME EIGENVALUES AND EIGENVECTORS OF T TO */
  /* SEE IF FURTHER ACTION IS INDICATED. */

  L310: if (*nop >=  *maxop || j >  *maxj -  *nblock)
  {
    test = TRUE_;
  }

  /* FIRST COMPUTE ALL THE EIGENVALUES LESS THAN AXL AND ONE MORE */

  L320: nl = numl;
  m = number + numl;
  dlaeig_(&j, nband, &c__1, &numl, &t[t_offset], &val[number + 1], maxj,
    &s[s_offset], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin, &tmax);
  if (val[m] > axl || numl == j)
  {
    goto L340;
  }
  ++numl;
  if (numl + numr >  *maxvec + 2)
  {
    goto L330;
  }

  i__1 = numr;
  for (ii = 1; ii <= i__1; ++ii)
  {
    i__ = numl + numr - ii;
    dcopy_(&j, &s[i__ *s_dim1 + 1], &c__1, &s[(i__ + 1) *s_dim1 + 1], &c__1);
    /* L325: */
  }
  dcopy_(&j, dzero, &c__0, &s[numl *s_dim1 + 1], &c__1);
  goto L320;

  /* FIND OUT HOW MANY THERE WHEN THERE ARE TOO MANY TO STORE */

  L330: nl =  *maxvec + 2-numr;
  m = number + nl;
  dcopy_(&j, dzero, &c__0, &s[nl *s_dim1 + 1], &c__1);
  dlaeig_(&j, nband, &numl, &numl, &t[t_offset], &val[m], maxj, &s[nl *s_dim1 +
    1], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin,  &tmax);
  if (val[m] > axl || numl == j)
  {
    goto L340;
  }
  ++numl;
  dcopy_(&j, dzero, &c__0, &s[nl *s_dim1 + 1], &c__1);
  goto L330;

  /* NOW COMPUTE ALL THE EIGENVALUES GREATER THAN AXR AND ONE MORE */

  L340: m = number + nl + 1;
  k = j - numr + 1;
  dlaeig_(&j, nband, &k, &j, &t[t_offset], &val[m], maxj, &s[(nl + 1) *s_dim1 +
    1], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin,  &tmax);
  if (val[m] < axr || numr == j)
  {
    goto L370;
  }
  ++numr;
  if (numl + numr >  *maxvec + 2)
  {
    goto L360;
  }
  i__1 = numr;
  for (ii = 2; ii <= i__1; ++ii)
  {
    i__ = numl + numr - ii + 1;
    dcopy_(&j, &s[i__ *s_dim1 + 1], &c__1, &s[(i__ + 1) *s_dim1 + 1], &c__1);
    /* L350: */
  }
  dcopy_(&j, dzero, &c__0, &s[(nl + 1) *s_dim1 + 1], &c__1);
  goto L340;

  /* FIND OUT HOW MANY THERE ARE WHEN THERE ARE TOO MANY TO STORE */

  L360: k = j - numr + 1;
  dlaeig_(&j, nband, &k, &k, &t[t_offset], &val[m], maxj, &s[(nl + 1) *s_dim1 +
    1], &bound[1], &atemp[1], &d__[1], &vtemp[1], eps, &tmin,  &tmax);
  if (val[m] < axr || numr == j)
  {
    goto L370;
  }
  dcopy_(&j, dzero, &c__0, &s[(nl + 1) *s_dim1 + 1], &c__1);
  ++numr;
  goto L360;

  /* THIS COMPUTES THE CORRESPONDING RESIDUAL NORMS AND */
  /* ORTHOGONALITY COEFFICIENTS. */

  L370: 
  /* Computing MIN */
  i__1 = numl + numr, i__2 =  *maxvec + 2;
  ntheta = min(i__1, i__2);
  dmvpc_(nblock, &bet[bet_offset], maxj, &j, &s[s_offset], &ntheta, &atemp[1],
    &vtemp[1], &d__[1]);

  /*  THIS CHECKS FOR TOO MANY VECTORS */

  if (numl + numr <=  *maxvec + 2)
  {
    goto L380;
  }
  *ierr = (numl + numr -  *maxvec - 2) * - 100;
  done = FALSE_;
  test = TRUE_;
  goto L440;

  /* THIS CHECKS IF TERMINATION IS ALLOWED */

  L380: m = number + nl;
  *deltal = val[m];
  *deltar = val[m + 1];
  if (j <  *nblock * 6 ||  *deltal - atemp[nl] < axl ||  *deltar + atemp[nl +
    1] > axr)
  {
    done = FALSE_;
  }
  if (done && ntheta == 2)
  {
    return 0;
  }
  if (ntheta == 2 && !test)
  {
    goto L80;
  }

  /* ------------------------------------------------------------------ */

  /* THIS SECTION EXAMINES THE COMPUTED EIGENPAIRS IN DETAIL. */

  /* THIS CHECKS WHETHER ALL THE EIGENVALUES ARE ACCEPTABLE. */

  L440: if (!(test || done))
  {
    goto L470;
  }
  nstart = 2;
  m = number + ntheta;
  /* Computing MAX */
  d__1 = rnorm, d__2 = val[m], d__1 = max(d__1, d__2), d__2 = val[number + 1];
  tola = utol * max(d__1, d__2);
  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (i__ == nl || i__ == nl + 1)
    {
      goto L460;
    }
    m = number + i__;
    /* Computing MIN */
    /* Computing MAX */
    d__3 =  *deltal - val[m], d__4 = val[m] -  *deltar;
    d__1 = atemp[i__] *atemp[i__] / max(d__3, d__4), d__2 = atemp[i__];
    if (min(d__1, d__2) > tola)
    {
      goto L450;
    }
    ind[i__] =  - 1;
    goto L460;

    L450: done = FALSE_;
    if (!test)
    {
      goto L470;
    }
    ++nstart;
    ind[i__] = 1;
    L460: ;
  }

  /* COPY IND VALUES INTO VTEMP */

  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (i__ == nl || i__ == nl + 1)
    {
      goto L465;
    }
    vtemp[i__] = (doublereal)((real)ind[i__]);
    L465: ;
  }
  goto L500;

  /* THIS CHECKS FOR NEW GOOD VECTORS. */

  L470: ng = 0;
  tolg = epsrt * anorm;
  i__1 = ntheta;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    if (i__ == nl || i__ == nl + 1)
    {
      goto L490;
    }
    if (vtemp[i__] > tolg)
    {
      goto L480;
    }
    ++ng;
    vtemp[i__] =  - 1.;
    goto L490;
    L480: vtemp[i__] = 1.;
    L490: ;
  }
  if (ng <= ngood)
  {
    goto L80;
  }
  nstart = ntheta - ng;

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES AND NORMALIZES THE INDICATED RITZ VECTORS. */
  /* IF NEEDED (TEST = .TRUE.), NEW STARTING VECTORS ARE COMPUTED. */

  L500: test = test && !done;
  /* Computing MIN */
  i__1 = ntheta - nstart, i__2 =  *maxvec -  *nperm;
  ngood = min(i__1, i__2);
  nstart = ntheta - ngood;
  vtemp[nl] = 1.;
  vtemp[nl + 1] = 1.;

  /* THIS ALIGNS THE DESIRED (ACCEPTABLE OR GOOD) EIGENVALUES AND */
  /* EIGENVECTORS OF T.  THE OTHER EIGENVECTORS ARE SAVED FOR */
  /* FORMING STARTING VECTORS, IF NECESSARY.  IT ALSO SHIFTS THE */
  /* EIGENVALUES TO OVERWRITE THE GOOD VALUES FROM THE PREVIOUS */
  /* PAUSE. */

  dcopy_(&ntheta, &val[number + 1], &c__1, &val[ *nperm + 1], &c__1);
  if (nstart == ntheta)
  {
    goto L530;
  }
  dvsort_(&ntheta, &vtemp[1], &atemp[1], &c__1, &val[ *nperm + 1], maxj, &j, 
    &s[s_offset]);

  /* THIS ACCUMULATES THE J-VECTORS USED TO FORM THE STARTING */
  /* VECTORS. */
  /* COMPUTE MINIMUM VALUE OF ATEMP TO AVOID POSSIBLE OVERFLOW */

  L530: if (!test)
  {
    nstart = 0;
  }
  if (!test)
  {
    goto L580;
  }
  temp = atemp[1];
  i__1 = nstart;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    /* Computing MIN */
    d__1 = temp, d__2 = atemp[i__];
    temp = min(d__1, d__2);
    /* L535: */
  }
  m = ngood + 1;
  l = ngood + min(nstart,  *nblock);
  i__1 = l;
  for (i__ = m; i__ <= i__1; ++i__)
  {
    d__1 = temp / atemp[i__];
    dscal_(&j, &d__1, &s[i__ *s_dim1 + 1], &c__1);
    /* L540: */
  }
  m = (nstart - 1) /  *nblock;
  if (m == 0)
  {
    goto L570;
  }
  l = ngood +  *nblock;
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nblock;
    for (k = 1; k <= i__2; ++k)
    {
      ++l;
      if (l > ntheta)
      {
        goto L570;
      }
      i1 = ngood + k;
      d__1 = temp / atemp[l];
      daxpy_(&j, &d__1, &s[l *s_dim1 + 1], &c__1, &s[i1 *s_dim1 + 1],  &c__1);
      /* L550: */
    }
    /* L560: */
  }
  L570: nstart = min(nstart,  *nblock);

  /* THIS STORES THE RESIDUAL NORMS OF THE NEW PERMANENT VECTORS. */

  L580: if (ngood == 0 || !(test || done))
  {
    goto L610;
  }
  i__1 = ngood;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    m =  *nperm + i__;
    res[m] = atemp[i__];
    /* L590: */
  }

  /* THIS COMPUTES THE RITZ VECTORS BY SEQUENTIALLY RECALLING THE */
  /* LANCZOS VECTORS. */

  L610: number =  *nperm + ngood;
  if (test || done)
  {
    i__1 =  *n **nblock;
    dcopy_(&i__1, dzero, &c__0, &p1[p1_offset], &c__1);
  }
  nstart = min(*nblock, nstart);
  if (ngood == 0)
  {
    goto L630;
  }
  m =  *nperm + 1;
  i__1 = number;
  for (i__ = m; i__ <= i__1; ++i__)
  {
    dcopy_(n, dzero, &c__0, &vec[i__ *vec_dim1 + 1], &c__1);
    /* L620: */
  }
  L630: i__1 = j;
  i__2 =  *nblock;
  for (i__ =  *nblock; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    (*iovect)(n, nblock, &p2[p2_offset], &i__, &c__1);
    i__3 =  *nblock;
    for (k = 1; k <= i__3; ++k)
    {
      m = i__ -  *nblock + k;
      if (nstart == 0)
      {
        goto L650;
      }
      i__4 = nstart;
      for (l = 1; l <= i__4; ++l)
      {
        i1 = ngood + l;
        daxpy_(n, &s[m + i1 * s_dim1], &p2[k *p2_dim1 + 1], &c__1, &p1[l
          *p1_dim1 + 1], &c__1);
        /* L640: */
      }
      L650: if (ngood == 0)
      {
        goto L670;
      }
      i__4 = ngood;
      for (l = 1; l <= i__4; ++l)
      {
        i1 = l +  *nperm;
        daxpy_(n, &s[m + l * s_dim1], &p2[k *p2_dim1 + 1], &c__1, &vec[i1
          *vec_dim1 + 1], &c__1);
        /* L660: */
      }
      L670: ;
    }
    /* L680: */
  }

  if (test || done)
  {
    goto L700;
  }

  /* THIS NORMALIZES THE RITZ VECTORS AND INITIALIZES THE */
  /* TAU RECURRENCE. */

  m =  *nperm + 1;
  i__2 = number;
  for (i__ = m; i__ <= i__2; ++i__)
  {
    temp = 1. / dnrm2_(n, &vec[i__ *vec_dim1 + 1], &c__1);
    dscal_(n, &temp, &vec[i__ *vec_dim1 + 1], &c__1);
    tau[i__] = 1.;
    otau[i__] = 1.;
    /* L690: */
  }
  goto L80;

  /* ------------------------------------------------------------------ */

  /* THIS SET IERR IF NECESSARY */

  L700: if (ntheta != 2 &&  *nop >=  *maxop)
  {
    *ierr +=  - 2;
  }
  if (ntheta == 2 &&  *nop >=  *maxop)
  {
    *ierr =  - 6;
  }
  if (ntheta == 2 && j >  *maxj -  *nblock &&  *ierr == 0)
  {
    *ierr =  - 5;
  }

  /* THIS SECTION PREPARES TO ITERATE THE ALGORITHM BY SORTING THE */
  /* PERMANENT VALUES, RESETTING SOME PARAMETERS, AND ORTHONORMALIZING */
  /* THE PERMANENT VECTORS. */

  if (ngood == 0 &&  *ierr == 0)
  {
    goto L30;
  }
  if (ngood == 0)
  {
    goto L800;
  }

  /* THIS ORTHONORMALIZES THE THE EIGENVECTORS */

  i__2 =  *nperm + ngood;
  dortqr_(nmvec, n, &i__2, &vec[vec_offset], &s[s_offset]);

  /* THIS SORTS THE VALUES AND VECTORS. */

  if (*nperm != 0)
  {
    i__2 =  *nperm + ngood;
    dvsort_(&i__2, &val[1], &res[1], &c__0, tarr, nmvec, n, &vec[vec_offset]);
  }
  *nperm += ngood;

  /* THIS DECIDES WHERE TO GO NEXT. */

  if (*ierr != 0)
  {
    goto L800;
  }
  if (!done)
  {
    goto L30;
  }

  /* THIS DOES A CLUSTER TEST TO SEE IF A CHECK RUN IS NEEDED */
  /* TO LOOK FOR UNDISCLOSED MULTIPLICITIES. */

  m =  *nperm -  *nblock + 1;
  if (m <= 0)
  {
    return 0;
  }
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    l = i__ +  *nblock - 1;
    if (val[l] - val[i__] < tola)
    {
      goto L30;
    }
    /* L790: */
  }

  /* THIS DOES A CLUSTER TEST TO SEE IF A FINAL RAYLEIGH-RITZ */
  /* PROCEDURE IS NEEDED. */

  L800: m =  *nperm -  *nblock;
  if (m <= 0)
  {
    return 0;
  }
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    l = i__ +  *nblock;
    if (val[l] - val[i__] >= tola)
    {
      goto L810;
    }
    *raritz = TRUE_;
    return 0;

    L810: ;
  }
  return 0;
} /* diwla_ */
